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Abstract 

The  linear  and  nonlinear  dynamic  behavior  of  flexibly-bladed  turbomachines  is 
presented.  The  analytical  description  is  based  on  matching  a  two  dimensional, 
incompressible  flow  field  across  a  semi-actuator  disk  representation  of  a  flexible 
rotor  and  a  rigid  stator.  The  aerodynamic  loading  on  the  rotor  is  derived  using 
control  volume  formulations  applied  to  discrete  blade  passages  allowing  consid¬ 
eration  of  finite  interblade  phase  angles.  Depending  on  operating  parameters, 
the  model  exhibits  behaviors  classified  as  surge,  rotating  stall,  and  stall  flutter 
which  are  qualitatively  consistent  with  experimentally  observed  results. 

The  formulation  provides  a  tractable,  nonlinear  state-space  description  of 
the  dynamics  responsible  for  surge,  rotating  stall,  flutter,  and  their  interaction. 
An  analysis  is  performed  for  system  parameters  representative  of  a  laboratory- 
scale  fan  test  facility.  The  behavior  of  the  operability-limiting  instability  modes 
is  examined  using  time  simulations,  eigenvalue  analysis  and  stability  diagrams. 
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1  Introduction 


Surge,  rotating  stall,  and  flutter  are  instability  phenomena  that  fundamentally 
constrain  the  design  and  operation  of  modern  gas  turbine  engines.  Classically, 
these  different  forms  of  compression  system  instabilities  have  been  categorized 
based  on  observations  of  observations  of  limit  cycle  oscillations.  Rotating  stall 
and  surge  are  primarily  aerodynamic  instabilities.  Flutter  is  an  aeroelastic  insta¬ 
bility  involving  blade  motion  which  can  occur  in  various  regions  of  the  compres¬ 
sor  performance  map.  In  addition  to  instability,  the  forced  response  of  flutter 
modes  is  extremely  important.  A  lightly-damped  flutter  mode  can  be  excited 
by  gusts  or  by  defects  in  the  inlet  flow  resulting  in  the  accumulation  of  fatigue. 

The  purpose  of  the  research  effort  is  to  provide  a  framework  to  develop 
and  analyze  nonlinear  dynamic  models  of  flexibly-bladed  turbomachines  where 
surge,  rotating  stall,  and  flutter  can  be  studied  as  coupled  instability  mecha¬ 
nisms.  One  cannot  expect  good  predictive  capability  from  such  a  model,  rather 
the  objective  is  to  determine  qualitatively  the  trends  in  system  stability  and 
nonlinear  dynamics  as  system  parameters  are  varied.  Reduced-order  models 
may  be  used  to  test  some  of  the  simplifications  typical  in  more  comprehensive 
computational  models  such  as  the  modeling  of  losses,  choice  of  the  inlet  bound¬ 
ary  condition,  and  the  influence  of  downstream  bladerows.  Another  role  for 
reduced-order  models  is  in  the  development  of  experimental  methodology  and 
in  the  interpretation  of  data.  Finally,  reduced-order  models  are  well  suited  for 
studies  in  active  and  passive  control. 

Reduced-order,  nonlinear  compression  system  models  have  been  constructed 
based  on  an  actuator  disk  approach.  A  hierarchy  of  models  has  been  developed 
which  includes  in  order  of  increasing  complexity:  (i)  a  compressor  model  with  an 
empirically-motivated  polynomial  characteristic,  (ii)  a  compressor  model  with 
a  characteristic  determined  from  the  ideal  pressure  rise  and  lagged  losses,  (iii) 
the  incorporation  of  an  inter-bladerow  space  between  the  rotor  and  stator,  and 
(iv)  the  incorporation  of  rotor  blade  flexibility.  The  models  are  in  the  form  of 
differential-algebraic  equations  constructed  from  idealized  governing  equations 
using  symbolic  algebra.  The  framework  is  such  that  changes  to  the  physical 
modeling  can  be  easily  incorporated. 

Analysis  has  focused  on  the  determination  of  stability  boundaries  through 
generalized  eigensystem  analyses  and  the  determination  of  the  initial  post¬ 
instability  dynamics.  These  two  features  provide  the  dominant  criteria  for 
predicting  the  operating  regime.  Principal  results  from  the  analysis  are  the 
demonstration  that  the  relative  damping  of  the  flutter  modes  is  determined 
most  strongly  by  the  plunge  to  twist  ratio  of  the  blades.  System  stability  has 
been  examined  with  respect  to  operating  conditions  including  shaft  speed,  nom¬ 
inal  mass  flow,  inlet  temperature  and  flight  speed.  It  is  also  shown  that  loss 
and  deviation  may  have  a  significant  effect  on  flutter  damping.  Numerical  sim¬ 
ulations  have  been  used  to  demonstrate  that  stable  flutter  and  stall  limit  cycles 
may  coexist.  Lastly,  linear  analysis  has  been  used  to  compute  pole-zero  dia- 
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grams  and  Bode  plots  demonstrating  how  one  would  use  the  models  to  support 
the  planning  of  experiments. 


2  Physical  modeling 

The  present  model  is  a  further  development  of  the  work  of  Gysling  and  My¬ 
ers  [2].  The  bladerows  are  represented  by  actuator  disks,  the  flow  fields  are 
two-dimensional  and  incompressible,  and  the  surge  dynamics  are  described  (ef¬ 
fectively)  by  a  Helmholtz  resonator.  The  system  variables  and  parameters  are 
listed  in  tables  1  and  2  respectively.  Wong  [8]  improved  upon  the  formulation 
of  Gysling  and  Myers  primarily  by  accounting  for  nonlinear  dependence  on  in¬ 
terblade  phase  angle.  The  present  model  differs  from  that  of  Gysling  and  Myers 
and  that  of  Wong  as  follows: 

•  In  the  formulation  of  the  deforming  rotor  control  volume,  it  is  shown  that 
the  width  of  the  control  volume  is  immaterial  as  one  should  expect  from 
an  actuator  disk.  In  those  earlier  works,  a  representative  blade  passage 
was  chosen  as  the  control  volume. 

•  The  present  model  accounts  for  control  volume  deformation  in  the  energy 
equation.  The  earlier  works  do  not. 

•  An  inter-bladerow  space  has  been  incorporated  for  the  purpose  of  param¬ 
eterizing  the  relative  coupling  between  the  rotor  and  the  stator. 

•  Inlet  total  pressure  is  introduced  as  a  parameter. 

•  Mass  flow  actuation  has  been  incorporated  for  the  purpose  of  supporting 
system  identification  and  control  experiments. 

The  flow  field  equations  are  linearized  about  a  condition  of  axisymmetric 
flow.  This  allows  us  to  solve  for  the  axial  dependence  of  the  non-axisymmetric 
flow  perturbation.  The  basic  variables  of  the  model  are  then  the  components  of 
velocity  and  pressure,  as  functions  of  the  azimuthal  angle  and  time,  at  five  dis¬ 
crete  axial  stations:  the  rotor  leading  edge,  the  rotor  trailing  edge,  just  upstream 
from  the  mass  flow  actuator,  just  downstream  from  the  mass  flow  actuator,  the 
stator  leading  edge,  and  the  stator  trailing  edge. 

Matching  conditions  across  the  actuator  disks  representing  the  bladerows 
are  given  by  the  general  expressions  for  the  conservation  of  mass,  momentum, 
and  energy.  The  equations  are  closed  using  strong  assumptions  about  flow 
kinematics,  i.e.  exit  flow  is  fixed  to  the  metal  angle  of  the  blade  plus  a  known 
deviation.  There  is  no  assumption  that  perturbations  are  small  and  the  resulting 
governing  equations  are  nonlinear  including  transcendental  functions  of  the  flow 
variables. 

The  key  assumptions  behind  the  model  are  incompressible,  two-dimensional 
modeling  of  the  flow,  actuator  disk  representation  of  the  bladerows,  and  the 
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Table  1:  List  of  system  variables 
<f>i  axial  velocity  at  the  rotor  leading  edge. 
v\  azimuthal  velocity  at  the  rotor  leading  edge. 
pi  pressure  at  the  rotor  leading  edge. 

<j}2a  axial  velocity  at  the  rotor  trailing  edge. 
v2 a  azimuthal  velocity  at  the  rotor  trailing  edge. 

p2a  pressure  at  the  rotor  trailing  edge. 

<j)2b  axial  velocity  just  upstream  of  the  mass  flow  actuator. 
v2b  azimuthal  velocity  just  upstream  of  the  mass  flow  actuator. 
p2b  pressure  just  upstream  of  the  mass  flow  actuator. 

(f>2c  axial  velocity  at  the  stator  leading  edge. 
v2c  azimuthal  velocity  at  the  stator  leading  edge. 
p2c  pressure  at  the  stator  leading  edge. 

<j>3  axial  velocity  at  the  stator  trailing  edge. 
pz  pressure  at  the  stator  trailing  edge. 
pp  pressure  in  the  exit  plenum. 

F  force  (vector)  exerted  on  the  blade  passage  by  the  flow. 

Lr  rotor  loss. 

Ls  stator  loss. 

q  blade  plunge. 

a  blade  twist. 

Bk  time  dependent  coefficient  in  inter-bladerow  space  streamfunction. 

Ck  time  dependent  coefficient  in  inter-bladerow  space  streamfunction. 

w  vorticity  at  rotor  trailing  edge. 

h  vorticity  lag. 


assumptions  on  flow  kinematics.  Wong  [8]  has  compared  predictions  of  a  similar 
model  to  two-dimensional  compressible  simulations  past  an  isolated  blade  row. 
He  found  that  reduced-order  model  was  suitable  for  stability  predictions  for 
Mach  number  less  than  0.8,  solidity  greater  than  8,  and  reduced  frequency  less 
than  0.15. 

2.1  Inlet  duct 

The  upstream  flow  field  is  assumed  to  be  irrotational  and  incompressible.  The 
system  is  assumed  to  draw  flow  from  a  reservoir  with  constant  total  pressure 
Patm.  Finite  inlet  effects  are  modeled  by  assuming  twodimensional  flow  in  a 
slab  whose  thickess  varies  axially  but  the  nominal  radius  of  which  is  constant 
as  suggested  by  Greitzer  [1],  A  sketch  of  the  idealized  geometry  and  the  flow 
control  volume  are  shown  in  figure  1.  The  governing  fluid  dynamical  equations 
are  the  conservation  of  mass  and  momentum. 

For  the  control  volume  sketched  in  figure  1,  conservation  of  mass  is  expressed 
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Table  2:  List  of  system  parameters 
atmospheric  pressure, 
inlet  geometry  parameter, 
actual  length  of  inlet  duct, 
actual  length  of  exit  duct, 
stator  chord, 
time  scale  for  rotor  loss, 
time  scale  for  stator  loss, 
trailing  edge  metal  angle  of  rotor, 
zero-incidence  angle  of  rotor  leading  edge, 
zero-incidence  angle  of  stator  leading  edge, 
coefficients  of  the  empirical  rotor  loss  function, 
coefficients  of  the  empirical  stator  loss  function.  . 
coefficients  of  the  empirical  rotor  deviation  function. 
Greitzer-B  parameter  (surge), 
plenum  loss  coefficient, 
throttle  parameter, 
length  of  inter-bladerow  space, 
number  of  blades 
rotor  chord. 

position  of  the  elastic  axis  divided  by  chord, 
position  of  the  center  of  gravity  divided  by  chord, 
position  of  the  center  of  pressure  divided  by  chord, 
rotational  inertia  divided  by  chord, 
blade  mass. 

frequency  of  pure  bending  (nondimension alized  by  ^shaft)* 

frequency  of  pure  twist  (nondimensionalized  by  ^shaft)- 

structural  damping  of  bending  mode. 

structural  damping  of  torsion  mode. 

stagger  angle  of  rotor. 

stagger  angle  of  stator. 

actuator  mass  flow. 
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Figure  1:  Sketch  of  idealized  inlet  and  control  volume. 


by  the  equation 

^-(PUh)Rnomdxd0+-^—^(pVh)Rnomdxd9  =  0, 

C/X  -rtnoin  Ou 

OI  Q  Q 

Txm+W(vh)  =  o, 

where  lengths  have  been  nondimensionalized  with  respect  to  the  nominal  radius. 
Considering  perturbations  about  an  axisymmetric  state, 

U  =  *(*,*) +  ^(M,0 

V  =  v(x,6,t) 

h  —  h(x), 


the  mass  conservation  equation  may  be  expressed  to  first  order 

d(j> ,  dh  dv 
teh  +  (t,te+hd9-0- 

For  an  irrotational  flow,  the  perturbations  may  be  expressed 

d_9F 

dx 
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V  = 


9F_ 

99’ 


and  the  mass  conservation  in  terms  of  the  velocity  potential  is  given  by 

,  _o  „  9F  dh 
h,V2F  +  -r-—  =  0. 
ox  ax 

If  h  is  constrained  to  depend  exponentially  on  x, 

h(x)  =  ~e~2aMttX, 

then  F  will  also  have  solutions  in  exponential  form, 

F  =  FneknXe™6  y 


where 

—  2ftinlet^n  “*  =  0. 

Clearly  the  perturbations  must  decay  as  x  — ►  — oo,  hence 

k„  =  QTinlet  +  yo^iet  +  n2‘ 

Applying  the  conditions  of  irrotationality  at  the  rotor  face,  <j> i  and  v\  are  related 
by  the  equations 

{kn  2c*inlet)<^ln  H"  i^ln  *=:  0.  (1) 

For  the  control  volume  sketched  in  figure  1,  conservation  of  momentum  is 
expressed  by  the  equation 

Q  rX+dx  /  1  \ 

—  J  pUhdx  -  -  \^p+ -pU2J  h, 

which  when  linearized  and  applied  at  the  rotor  face  yields 


<f>\n 


kn  2Qfiniet 


=  “(Pin  + 


(2) 


The  axisymmetric  variables  are  assumed  to  be  related  through  a  momentum 
equation  through  a  finite-length,  uniform  duct, 


(P* tm  -  j*?)  -  Pi  =  iinletil. 


(3) 
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2.2  Blade  passage  deformation 

The  blade  loading  and  matching  conditions  across  the  actuator  disk  are  derived 
by  considering  a  control  volume  corresponding  to  the  passage  between  two  ad¬ 
jacent  blades.  The  control  volume  will  deform  of  course  with  blade  deflection. 
Deformation  of  the  leading  and  trailing  edge  boundaries  directly  affects  the 
conservation  of  mass,  momentum,  and  energy  in  the  control  volume.  The  con¬ 
trol  volume  is  sketched  in  figure  2  and  the  positive  sense  of  blade  deflection  is 
indicated  in  figure  3. 


Figure  2:  Definition  sketch  for  blade  passage  control  volume. 

The  deflections  of  the  actuator  disk  are  defined  such  that  the  actual  blade 
deflections  may  be  obtained  simply  by  collocation.  The  twist  of  blade  i,  for 
example  is  given  by 

ot 

Choosing  the  undeflected  elastic  axis  of  the  first  blade  as  the  origin,  the 
coordinates  of  the  leading  and  trailing  edges  are  given  by 

X\e  =  -q  sin  jT  -  £eac  cos(yr  -  a) 

6\e  —  6  +  q  cos  7r  -  £eac  sin(7r  -  a) 

Xte  =  -?sin7r  +  (l-£ea)ccos(7r-a) 

6te  =  0-bgcos7r  +  (1 -£ea)csin(7r -- a). 
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Figure  3:  Definition  sketch  for  blade  deflection  indicating  positive  sense  of  twist 
a  and  plunge  q.  The  location  of  the  elastic  axis  is  indicated  by  the  circle. 


The  outward  normal  vectors  are  given  by 

nie  =  -  cos  +  sin  /?i«j 

nte  =  cos  Ptei  -  sin  /?t  J , 


where  the  angles  /?|e  and  /?te  are  given  by 


tan  p\e 
tan  pte 


dxle/dO 

eele/de 

dxte/d6 

de^/oe  ‘ 


The  path  lengths  along  the  leading  and  trailing  edges  are  given  by 


It  is  assumed  that  the  exit  angle  of  the  flow  is  fixed  by  the  metal  angle  of 
the  deflected  blade  plus  some  deviation.  The  deviation  and  also  the  quasi-static 


total  pressure  loss  are  assumed  to  be  functions  of  the  incidence  angle.  The  flow 
geometry  is  sketched  in  figure  4.  The  relative  velocity  between  the  flow  and  the 
leading  edge  is  given  by 

Vr,le  =  (<£li  +  t>lj)  ~  (Jjj  ~  -fijjj  (Zle*  +  ^lej)i 

The  rotor  incidence  angle  is  then 

Gjnc.r  =  tan"1  ( ^-4)  -  /?zr  + 

Vvr,le*l/ 

where  /?zr  is  the  zero  incidence  angle  of  the  undeflected  blade.  Similarly,  the 
relative  velocity  between  the  flow  and  the  trailing  edge  is 

(d  d  \ 
dt  ~  86  J  + 

The  kinematic  constraint  on  the  exit  flow  angle  is  then 


-  (l  -  w) e “ = h*  -  (I  - 1)  *-] tanw'  -  ° + a)i 


(4) 


where  /?r  is  the  metal  angle  of  the  trailing  edge  and  the  deviation  is  assumed  to 
be  given  by 

A  =  AiOfinc, r  *4*  A2. 


2.3  Conservation  of  mass  across  the  rotor 

Conservation  of  mass  within  the  deforming  control  volume  may  be  expressed  as 
follows: 


^  J  pdV  +  J p(vT-n)dS  =  0. 


The  subscript  r  on  the  velocity  vector  indicates  velocity  with  respect  to  the 
moving  boundary.  The  time-derivative  on  the  volume  integral  is  to  be  taken 
in  the  rotor  frame.  The  present  model  is  two-dimensional,  incompressible,  and 
mass  flux  occurs  over  the  leading  and  trailing-edge  boundaries  only,  hence  the 
mass  conservation  equation  may  be  expressed 


' 8 
y8t 

or 

a 

8  > 
86) 

\  9V  (  v  <9sie  ,  \  <9«Ste 

1  -gj  +  (vr,le-nie)  -jp  +  (vr,te-nte)  QQ 

where 

8V 

86 

(  Jfds\e  8ste\ 

=  -  “>2  +  m  ■ 

=  0, 


(5) 


13 


vr,le 


Figure  4:  Geometry  of  leading  and  trailing-edge  flow. 


2.4  Conservation  of  momentum  across  the  rotor 

Conservation  of  momentum  within  a  control  volume  may  be  expressed 

—  j  pvdV  +  J px  (vr  •  n)  dS  -f  J pndS  =  — F. 

Applying  this  equation  to  the  deforming  blade  passage  yields 

(I  - §>)  sy> + £ (v  (v'  ■ n) + H  %ds =-£  f 1 "■ 


or 


(I  -  Te)  (vw)  +  [V|* (Vr'" 'n") + 

'  r  ,  v  ,  ,  dste  dF 

+  [vte  (vr,te-nte)  +  Ptente]  -fig-  =  —fig, 


de 


(6) 


where 


=/ 

Je, 


9 3  dF 

le 


de 
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is  the  force  exerted  on  the  control  volume  boundaries  by  the  fluid  and 

vie  =  ^ii  +  (l  +  Vi)j 

Vte  =  </>2ai+  (1  +  Via)} 

Pie  =  Pi 
Pte  P2a- 

Additionally,  the  velocity  within  the  control  volume  has  been  approximated  by 
the  average  of  the  leading  and  trailing-edge  flow  velocities, 

v  =  i(vle  +  Vte). 


2.5  Conservation  of  energy  across  the  rotor 


Conservation  of  energy  within  a  control  volume  may  be  expressed 

lJ1-Pv*dV  +  j  (p+ipv2)(vr.n)dS  =  -F.v  +  C?. 
Applying  this  equation  to  the  deforming  blade  passage  yields 

- 1)  f  + i,”  (* + r?-  - Lr )  (v'J'ni*) 

( d  d\  1  2dv  ,  /  ,1  2  r  \  ,  ,ds]e 

(di-Te)  2V  W  +  (Ple  +  2^  "  Lr)  (Vr'le'nie)  ~ef 

+  ^Pte  +  2V‘e^  (Vr.te-nte)  =  ~~QQ  '  Vcv> 


where  Lr  represents  a  loss  in  leading-edge  total  pressure  to  account  for  noncon¬ 
servative  processes.  The  velocity  of  the  control  volume  is  approximated  by  the 
average  of  the  velocities  of  the  leading  and  trailing-edge  boundaries. 


v"  =  (i 


d  \  ( X[e  +  Xte  .  01e  +  #te  A  . 

de)\  2  2  v  J 


2.6  Blade  dynamics 

The  blade  dynamics  are  modeled  using  a  typical  section  with  inertial  and  aero¬ 
dynamic  coupling  between  twist  and  plunge.  The  lift  force  is  assumed  to  act  at 
the  center  of  pressure  which  is  assumed  constant, 

Fi  =  [F(0  +  n/Nb,t)  -  F(0  -  7r /Nb,t)]  .[-sin(7r  -  <*)i  +  cos(7r  -  a)j], 
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and  the  aerodynamic  moment  is  applied  at  the  elastic  axis, 
constant, 


M  —  Fl(£e a  ”  £cp)c. 


which  is  also  assumed 


The  blade  bending  equation  is  then 

+2CbQb  -  gfij  q  +  Qti  =  p.  (8) 


and  the  blade  twist  equation  is 


where 


&  4-  (^ea  ~  &g)c-P 
dO)  +  Jea 


+2CtQt  (jt 


/ea  =  D(U  -  6g)2c2  +  £>cV. 


(9) 


2.7  Flow  field  in  the  inter-bladerow  space 

In  the  inter-blade  space,  the  perturbation  to  the  axisymmetric  flow  is  modeled 
as  a  vortical  inviscid  flow.  The  axisymmetric  components  of  velocity  V  and  $ 
remain  constant  over  the  interblade  space, 

$2a  =  $26  (10) 

v2a  =  Vu;  (11) 


the  axisymmetric  component  of  the  pressure  field  may  be  found  from  the  balance 
of  energy, 

P2a  -  P26  =  libs  (12) 

The  velocity  and  pressure  perturbation  fields  at  the  upstream  and  down¬ 
stream  ends  of  the  inter-bladerow  space  will  be  related  to  each  other  through 
the  streamfunction  #  and  the  associated  vorticity  field  w, 


-V2*  = 
Du  __ 
~Dt 


w, 

du ;  _  du>  T_  du>  n 

!h+®2a!h  +  V2a7)e  _0' 


The  perturbation  streamfunction  can  be  decomposed  into  homogenous  and 
particular  components  -1-  \Pp,  such  that 


—  0,  and  —  V2\Pp  =  u. 
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The  general  solution  for  the  homogenous  component  can  easily  be  shown  to  be 

$h(x,  9,  t )  =  Bk(t)  exp  (kx  +  ik9)  +  Ct(t)  exp(-fcx  +  ikO). 

An  approximation  for  the  particular  solution  follows. 

The  vorticity  convection  equation  has  the  solution 


If  harmonic  dependence  on  9  is  required  and  if  u>o,k(t)  is  defined  such  that 

wCO,#,*)  =  w0}k(t)ex?(ikO), 

"<*•  = "»■*  ('  -  sb) exp  [*  (*  -  Iff)]  - 

The  convective  delay  will  be  approximated  by  a  lag  equation.  If  u>o  and  ujx 
are  defined  such  that 

w0(t)  —  0)  and  ujx(t^  —  6 ,  t), 

then  clearly  wx(t  +  r)  =  wo(0>  where  T  —  x/$. 

The  Laplace  transform  of  uo  is  given  by 

C[u)o(t)]  =  C[u>x(t  +  t)]  =  eSTC[u>L(t)  ]. 

Hence, 

^i^x)  _  -5T 

C(u o)  "  * 

If  the  exponential  is  replaced  by  its  (1-1)  Pade  approximation,  one  obtains 
C(ux) 

C(u)  o) 

Introducing  h(t)  such  that 


yields  the  two  equations: 

(2  +  sr)C(h)  =  4£(u/0), 

C(ujx)  =  — £(wo)  +  C{K). 

Applying  the  inverse  Laplace  transformation,  one  may  solve  for  u>x  in  terms  of 
u>o  and  h, 

wx(t)  =  h(t)-wo(t) 

dh  2  ' .  » 

_  =  _(_A+2wo). 
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=  -i  +  — - — . 
2  +  sr  2  +  sr 


C{h)  = 


4£(t>>o) 
2  4-  st  ’ 


Thus  for  x  >  0, 

u{x ,  6,  t)  «  [hk(x,  t)  -  w0>*(t)]  exp (ik0  -  ikV */$), 

(13) 

ot  X 

Finally,  the  dependence  of  the  vorticity  on  x  must  be  considered.  If  the 
vorticity  is  assumed  to  vary  linearly  with  x, 


u >(x,9,t)  =  ^1  -  — ^  w(0 .dyt)  + 


V  17  '  " 

=  (l  -  |)  u>o,t(t)exp(ik0)  +  j[hk(t )  -  u0,k(t)]exp{ik9  -  ikVL/9), 


where 


dhjs  2$>  ,  0 

_  =  T(-h*  +  2u,0l*), 


then  the  particular  component  of  the  streamfunction  must  also  vary  linearly 
with  x,  hence, 

p{x,e,t )  =  j  J u{x,6,t)d6dd, 

=  p-  {  (l  -  ^ )  wo,k(0  exp(iM) 

+  j[hk(t)  -wo,fc(0]exp(iA:ff  +  t*VI/<&)}  • 

The  boundary  conditions  at  the  upstream  end  of  the  inter-blade  space  to¬ 
gether  with  equation  13  relate  the  constants  defining  the  streamfunction,  Bk(t), 
Ck(t),  u>o, k{t),  and  hk{t)  to  the  perturbation  quantities  <t>2a(0,i ),  V2a(9,t),P2a(9,t) 
and  the  nominal  velocity  components  <f>2a  and  Vja : 

<t>2  a{9,t)  =  +  -QQ  ^  (^) 

V2a(9,t)  -  ~  -fa  (15) 

_  \*±  +  **l  (16) 

de  ~  dxdt +  dx 2  dxde J  x=0 

Likewise,  the  boundary  conditions  at  the  downstream  end  relate  the  perturba¬ 
tion  quantities  V2i(0,O.P2»(M)  to  the  constants  defining  the  stream- 

function  Bjfc(t),  C'fc(t),  wq.j fc(t),  hk(t),  and  $2a  and  V2a> 


4>2b{9,t) 


*  U. 
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V2b(0,t) 


(18) 


aw 


dp2b 

de 

2.8  Actuation  by  bleed  of  mass  flow 

The  bleeding  of  mass  flow  is  a  physically-realizable  means  of  control.  Bleed  is 
incorporated  into  the  present  reduced-order  model  through  an  actuator  disk  just 
downstream  from  the  inter-bladerow  space.  As  one  would  expect,  the  efficacy 
of  control  is  strongly  dependent  on  the  axial  position  of  the  bleed  valve. 

The  actuator  disk  matching  conditions  are  found  from  the  balance  of  mass 
and  momentum  flow  through  the  control  volume  sketched  in  figure  (5).  The 


H”  ^  Q  O  V  C\ 


^blced 


Figure  5:  Sketch  of  bleed  control  volume. 


axial  dimension  of  the  control  volume  is  taken  to  be  sufficiently  small  that  the 
instantaneous  mass  and  momentum  within  may  be  neglected.  The  mass  flow 
equation  is 

fhjn  rhbieed  ^out  ~  0> 
the  equation  for  the  flow  of  axial  momentum  is 


-^in^  P 


dt(mu) 
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^in^in  ^bleed^bleed  ^out^out) 


and  the  equation  for  the  flow  of  tangential  momentum  is  likewise 

0  =  |(rn») 

=  ^in  ^in  ^^bleed  ^bleed  ^hout^out* 

As  an  additional  constraint,  it  is  assumed  that  Ubieed  =  Ujn  and  Vbieed  =  ^in- 
The  axial  dimension  of  the  control  volume  is  assumed  to  sufficiently  small  that 
the  instantaneous  mass  within  can  be  neglected. 

The  non  dimensional  equations  for  the  bleed  actuator  disk  are  found  to  be 


02*  —  £  “  02c  =  0 

(20) 

P2b  ~  P2c  -  £026  4“  £2  =  0 

(21) 

V2 b  -  V2c  =  0, 

(22) 

where  all  quantities  depend  on  9  and  t. 


2.9  Conservation  equations  for  the  stator 

The  stator  is  modeled  as  a  rigid  blade  row.  Matching  conditions  for  the  actuator 
disk  are  given  by  conservation  of  mass,  the  Bernoulli  equation  with  the  inclusion 
of  a  loss  term.  Further  it  is  assumed  that  the  exit  angle  of  the  flow  is  fixed  by 
the  metal  angle  of  the  blade  (which  is  taken  to  be  zero).  Conservation  of  mass 
can  be  expressed 

02c  =  <t>  35  (23) 

the  Bernoulli  equation  for  the  flow  across  the  stator  is 


cs  d(j)2c 

cos  7S  dt 


(24) 


where  Ls  represents  a  loss  in  leading-edge  total  pressure  and  V3  is  fixed  to  be 
zero. 


2.10  Lagged-loss  equations 

The  total  pressure  losses  across  the  rotor  and  stator  actuator  disks  are  assumed 
to  lag  their  quasi-static  values.  A  simple  one-dimensional  lag  equation  is  used 
in  each  case, 

Tr  (J? "  §b)  Lt  =  ~iLr  ~  Ir'qs)’  (25) 

and 

n  T 

(26) 


dL»  -  (T  T  '> 

7s  ~ .  —  \Ls  L/Siqs/' 


dt 


20 


The  quasi-static  losses  are  taken  to  be  functions  of  incidence  angle, 

£r,qs  =  -krlttfncjr  +  -^r2^inc,r  +  ^r3 
Ls,qs  =  ^sl^fnCjS  +  ^s2ainc,s  +  ^s3> 

where 

ainc.s  =  -  tan-1  “  A»- 

Ofinc.r  is  defined  in  section  2.2. 


2.11  Exit  duct 

The  flow  in  the  exit  duct  is  assumed  to  be  incompressible  and  rotational.  It  is 
assumed  that  the  stator  fixes  the  exit  flow  angle  to  zero  and  the  length  of  the 
duct  is  assumed  infinite.  The  more  general  case  of  nonzero  exit  flow  angle  is 
treated  in  reference  [4]. 

To  first  order,  the  pressure  perturbation  satisfies  a  Poisson  equation, 


V2p  =  o, 


which  yields  solutions  of  the  form 


—  nx  in 9 


P  =  PSn(t)e  nxe] 

The  assumption  of  an  infinite  duct  precludes  the  possibility  of  upstream  prop¬ 
agating  waves.  The  flow  must  also  satisfy  the  linearized  axial  momentum  and 
continuity  equations, 

d<fr  ^  9(j)  _ 

8t  dx^  80  F’ 


and 


dv 

8x  +  80 


At  the  exit  face  of  the  stator  (*  =  0)  the  flow  angle  is  assumed  to  be  zero, 
thus, 


Continuity  requires  then 


d<j> 

dx 


=  0, 


at  x  =  0,  thus  the  axial  momentum  equation  reduces  to 

<P  3n  =  «P3n- 


(27) 


As  in  the  inlet,  the  axisymmetric  variables  are  assumed  to  be  related  through 
a  momentum  equation  through  a  finite-length,  uniform  duct.  Following  Moore 


21 


and  Greitzer  [6],  the  exit  duct  is  assumed  to  empty  into  a  plenum.  From  the 
plenum,  flow  exits  through  an  orifice  back  into  the  constant  total  pressure  reser¬ 
voir.  Moore  and  Greitzer  were  motivated  by  the  stability  of  compressors  for 
which  it  is  reasonable  to  assume  a  complete  loss  of  dynamic  head  when  the  flow 
dumps  into  the  plenum.  The  present  analysis  is  motivated  by  the  bypass  flow 
of  a  large  turbofan  engine.  In  this  case  the  exit  duct  simply  discharges  into  the 
atmosphere. 

The  equation  for  momentum  conservation  in  the  duct  is  taken  to  be 

(ft  +  5*1)  -  (ft  +  5 CM)  =  iexit*3.  (28) 

As  a  generalization  of  the  surge  modeling  of  Moore  and  Greitzer,  Cd  is  intro¬ 
duced  as  a  coefficient  modifying  the  loss  of  dynamic  head.  The  extreme  cases 
are  given  by  Cd  =  1  in  which  all  dynamic  head  is  lost  as  the  duct  dumps  into  the 
plenum  and  by  Cd  =  0  in  which  the  flow  passes  without  loss  through  a  diffuser. 

The  equation  for  mass  conservation  within  the  plenum  may  be  expressed 

4 B2LePp  =  *3  -  Ay/tyPp  -  fttm),  (29) 

following  Moore  and  Greitzer.  Greitzer’s  B  parameter  is  proportional  to  the 
square  root  of  the  volume  of  the  plenum.  The  bypass  flow  of  a  turbofan  is  best 
modeled  by  the  combination 

Cd  =  B  =  0. 
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3  Model  reduction 


Each  system  variable  is  projected  over  a  finite  number  of  Fourier  components 
in  the  azimuthal  variable.  For  example, 

N 

<j>i(0,t)  =  <j>i,i{t)  +  <£i,2n  cos 710  -  <j>i,2n+i  s\nn6. 

n=l 

Thus,  in  equation  3, 
and  in  1 

<f)ln(t)  -  +  tyl,2n+l(t)' 

A  time-dependent  system  of  equations  is  derived  through  a  Galerkin  projec¬ 
tion  of  equations  1  —  29  over  the  Fourier  modes.  Equations  3,10,11,12,28,  and 
29  relate  the  nominal  axisymmetric  flow  quantities  and  are  already  functions  of 
time  alone.  Equations  1,  2,  13,  14,  15,  16,  17,  18,  19,  and  27  relate  the  flow 
perturbations  only  and  need  must  be  projected  over  the  basis  {cos  n0,  sin  n9), 
for  n  =  1 . .  .JV.  Finally,  equations  4,  5,  6,  7,  8,  9,  20,  21,  22,  23,  24,  25,  and  26 
involve  axisymmetric  and  perturbation  terms  and  thus  must  be  projected  over 
the  basis  {1,  cos  nO ,  sin  n6 },  for  n  —  1 . . .  N.  Equations  8  and  9  are  second-order 
in  time,  and  equation  6  is  a  vector  equation  with  two  components.  Written  as 
a  system  of  first  order  differential-algebraic  equations, 

/(*,£, p)  =  0, 

the  total  number  of  equations  and  state  variables  is  52/S7  -f  22. 
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4  Linear  stability  analysis 

Linear  stability  is  a  key  element  in  determining  the  bounds  on  the  parameter 
space  for  design  and  operation  of  the  turbomachine.  It  is  feasible  using  the 
reduced-order  model  to  quickly  determine  system  damping  and  qualitatively- 
accurate  boundaries  to  linear  stability  as  the  system  parameters  are  varied.  This 
will  be  the  primary  tool  for  predicting  the  system  operability  map,  determining 
the  sensitivity  of  flutter  instability  to  installation  effects,  viscosity  (losses  and 
deviation),  and  mistuning,  and  investigating  system  identification  and  control. 

The  reduced-order  model  is  expressed  as  an  autonomous  system  of  differential- 
algebraic  equations 

A(«,i,p)  =  0,  (30) 

where  x  is  the  state  vector  and  p  is  a  vector  of  parameters.  The  fixed  points  xp 
of  system  (30)  satisfy  the  nonlinear  algebraic  system 

A(x,0,p)  =  0, 

which  may  be  solved  iteratively.  The  locus  of  fixed  points  as  the  throttle  is 
varied  at  constant  shaft  speed  correspond  to  a  characteristic  or  speedline  of  the 
turbomachine.  Typically  plotted  as  pressure  rise  versus  mass  flow,  the  charac¬ 
teristic  has  a  peak  at  which  the  maximum  possible  pressure  rise  is  achieved. 

The  stability  of  each  point  on  the  speedline  may  be  determined  from  the 
generalized  eigenvalue  problem 

[XAx  4-  Ax]x  =  0,  (31) 


where 


and  Ax 


the  linearization  of  system  (30)  about  the  fixed  point  x.  The  eigenvalue  problem 
(31)  is  solved  numerically  using  the  RGG  routine  of  the  Slatec  (Eispack) 
subroutine  library  [7].  The  eigenvalues  A  are  returned  as  a  ratio 


A  =  a//?, 


where  a  is  complex  and  /?  is  real.  The  matrix  A*  is  typically  singular,  hence, 
only  the  eigenvalues  for  which  (3  is  nonzero  are  relevant  to  system  stability. 

The  relative  stability  of  a  system  eigenmode  is  concisely  described  by  its 
damping  ratio.  If  the  eigenvalue  is  given  by 


the  damping  ratio  is 


A  =  a  4-  iu>, 


<  =  - 


a 

y/u)2  +  a2 
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(  <  0  indicates  the  mode  is  unstable.  The  damping  of  flutter  modes  is  computed 
with  respect  to  the  blade-fixed  coordinates,  hence 

Q 

yj  (w2  —  in)  +  or2  ’ 

where  n  is  the  Fourier  nodal  diameter  of  the  flutter  mode  (positive  or  negative). 


25 


5  System  stability 

The  parameters  governing  operation  are  the  throttle  parameter,  the  shaft  speed 
(by  which  time  is  nondimensionalized),  and  the  inlet  total  pressure.  At  constant 
altitude  operation,  such  as  in  a  test  facility,  it  is  reasonable  to  fix  the  inlet  total 
pressure  to  a  constant  value.  The  locus  of  equilibria  as  the  throttle  is  varied  at 
constant  shaft  speed  is  referred  to  as  a  speedline.  It  is  often  plotted  as  total- 
to-static  pressure  rise  against  the  flow  coefficient.  The  operating  line  or  “op 
line”  is  the  locus  of  equilibria  as  the  shaft  speed  is  varied  at  a  particular  value 
of  the  throttle  parameter.  In  this  section,  the  dependence  of  system  stability 
on  throttle  and  shaft  speed  is  considered.  All  other  parameters  are  assigned 
nominal  values  consistent  with  a  small-scale  research  rig.  These  are  given  in 
table  3.  The  dependence  of  system  stability  on  blade  geometry,  duct  lengths, 
and  other  parameters  is  considered  elsewhere. 

The  fundamental  system  instabilities  are  high-incidence  flutter  and  rotating 
stall.  These  occur  as  circumferentially  travelling  waves.  Flutter  can  be  either 
forward  or  backward  travelling.  Stall  is  always  a  forward  travelling  wave.  Due 
to  axial  symmetry,  the  instability  modes  decouple  into  their  Fourier  harmonics, 
or  nodal  diameter  patterns.  For  large  values  of  the  Greitzer-5  parameter  surge 
is  also  a  fundamental  instability.  The  flutter  eigenvalues  lie  very  near  their 
natural  frequencies  i(u;&  ±  n)  and  i(c ot  ±  n)  consistent  with  the  understanding 
that  the  aeroloading  is  small  in  comparison  to  the  elastic  restoring  forces. 

Neglecting  structural  damping  and  aero-loading  in  the  (rotor  frame)  blade 
deflection  equations  8  and  9,  yields  the  following  conservative  system: 

Q  4"  (£ea  ^cg)^  +  QbQ  =  0, 


&+iU'^)cDq  +  Qla  =  0. 

^ea 

Except  for  the  case  in  which  the  elastic  axis  and  center  of  gravity  are  coincident, 
each  of  the  two  mode  shapes  involves  coupled  bending  and  twist.  Assuming 
harmonic  response,  it  is  easily  shown  that  the  natural  frequencies  of  the  flutter 
modes  are  given  by 


U) 


2 

b 


Ut 


qi  +  qI-  VM  ~  W  ±  WM 

2(1  -  A*) 

Ql  +  Ql  +  ViQl-Qir  +  MiQT) 

2(1  -a0 


(32) 


where 


—  (£ea  ~  £°g)2 
(£ea  -  £cg)2  +  e2 


and  it  is  assumed  that  u>t  <  W(.  In  the  limit  n  — *  0,  w;,  =  Qt,  and  =  Qt . 
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Table  3:  Nominal  values  for  system  parameters. 


-^atm 

=  2.17  x  106  s~2/Qghaft 

atmospheric  pressure. 

cs/cos(7s) 

=  0.2086, 

effective  length  of  stator  blade  passage. 

fir 

=  15°, 

trailing  edge  metal  angle  of  rotor, 

Pzr 

=  45°, 

zero-incidence  angle  of  rotor  leading  edge, 

/?zs 

=  32°, 

zero-incidence  angle  of  stator  leading  edge, 

B 

=  0, 

Greitzer-5  parameter  (surge). 

Cd 

=  0, 

plenum  loss  coefficient. 

Nb 

=  16, 

number  of  blades 

cr 

=  0.4560, 

rotor  chord. 

D 

=  14.6, 

blade  mass. 

1 r 

II 

CO 

o 

0 

stagger  angle. 

£*inlet 

=  0.5, 

inlet  geometry  parameter. 

-^inlet 

=  0.5, 

actual  length  of  inlet  duct. 

^exit 

=  0.5, 

actual  length  of  exit  duct. 

Tr 

=  0.61, 

time  scale  for  rotor  loss. 

rs 

=  0.32, 

time  scale  for  stator  loss. 

Ai 

=  0.18, 

coefficient  of  the  empirical  rotor  deviation  function. 

a2 

=  12°, 

coefficient  of  the  empirical  rotor  deviation  function. 

■^ibs 

=  0.3, 

length  of  inter-bladerow  space.  ' 

fcg 

=  0.35, 

position  of  the  center  of  gravity  divided  by  chord. 

£cp 

=  0.35, 

position  of  the  center  of  pressure  divided  by  chord. 

Lri 

=  1.8842, 

coefficients  of  the  empirical  rotor  loss  function. 

Lt2 

=  -0.5053, 

LT3 

=  0.1219, 

Lsi 

=  0.7429, 

coefficients  of  the  empirical  stator  loss  function. 

Ls2 

=  -0.1450, 

Lsz 

=  0.0951. 

£ea 

=  0.55, 

position  of  the  elastic  axis  divided  by  chord. 

e 

=  0.20, 

rotational  inertia  divided  by  chord. 

Qb 

=  2620  s  ^/fishaft, 

frequency  of  pure  bending  (nondimensionalized  by  Qshaft)* 

Qt 

=  1990s-Vnshaft, 

frequency  of  pure  twist  (nondimensionalized  by  fishaft)* 

Cb 

=  0.035, 

structural  damping  of  bending  mode. 

Ci 

=  0.035, 

structural  damping  of  torsion  mode. 

27 


Solution  of  the  fixed  point  problem  and  generalized  eigensystem  analysis 
are  sufficient  to  determine  system  stability  but  not  system  operability  because 
rotating  stall  typically  occurs  as  a  subcritical  Hopf  bifurcation,  hence  there  is  a 
region  of  parameter  space  in  which  the  system  equilibrium  is  locally  stable  but 
globally  unstable. 

In  figure  (6)  is  plotted  a  speedline  for  fishaft  =  942s"1  (9000  RPM).  The 
fixed  point  value  of  plenum  pressure  is  plotted  against  the  fixed  point  value 
of  flow  coefficient  the  throttle  parameter  is  varied.  The  fixed  points  of  system 
30  are  computed  using  Newton- Raphson  iteration.  The  points  of  flutter  and 
stall  instability  are  indicated.  Stall  instability  typically  occurs  near  the  peak  of 
the  total-to-static  pressure  rise.  If  it  is  assumed  that  the  compression  system 
operates  in  a  quasi-steady  manner  such  that  the  pressure  rise  is  a  function  of 
flow  coefficient  only,  it  can  be  shown  that  instability  always  occurs  at  the  peak. 
In  reality,  compression  systems  have  a  finite  response  time  and  it  has  been  shown 
that  this  has  a  stabilizing  effect  [3]. 

Loci  of  stator-frame  eigenvalues  corresponding  to  the  speedline  of  figure  (6) 
are  plotted  in  figure  (7).  For  clarity,  only  a  subset  of  the  least  stable  modes 
is  shown.  The  flutter  and  stall  frequencies  are  nearly  constant  over  the  speed¬ 
line.  The  growth  rates  (determined  by  real  part)  are  very  different  however; 
in  comparison  to  the  flutter  eigenvalues,  the  stall  eigenvalue  moves  much  faster 
from  very  stable  to  instability  over  the  speedline.  The  variation  of  blade-frame 
damping  ratio  along  the  speedline  for  the  lowest  nodal  diameter  mode-1  flutter 
modes  is  plotted  in  figure  (8).  At  a  fixed  point  on  the  speedline,  the  damping  of 
the  mode  1  flutter  modes  is  plotted  in  figure  (9).  The  current  configuration  has 
eight  blades,  hence  the  x-axis  of  figure  (9)  corresponds  to  a  range  of  inter-blade 
phase  angle  from  -180  to  180  degrees.  The  variation  of  stall  damping  ratios 
over  the  speedline  is  plotted  in  figure  (10).  Note  that  the  scale  over  which  stall 
damping  varies  is  much  larger  than  the  scale  over  which  flutter  damping  varies. 
For  the  chosen  set  of  parameters,  the  initial  instability  is  forward-travelling  first 
nodal  diameter  flutter. 

The  operability  of  a  engine  can  be  characterized  by  its  response  as  the  mass 
flow  rate  and  shaft  speed  are  independently  varied.  Other  important  effects 
which  are  not  considered  here  are  nonuniform  inlet  flows,  including  flows  with 
nonzero  incidence  and  gusts.  The  sensitivity  to  inlet  total  pressure  is  considered 
in  a  later  section  (6.7).  The  stability  boundaries  of  the  least-stable  modes  as 
the  throttle  parameter  and  shaft  speed  are  independently  varied  are  plotted  in 
figure  11.  The  region  of  stable  operation  is  bounded  by  the  zero-nodal  and  first 
forward-travelling  flutter  boundaries.  For  this  set  of  parameters,  flutter  always 
occurs  before  stall.  Linear  stability  however  is  not  a  conservative  measure  of 
operability;  the  stall  instability  is  typically  subcritical,  thus  there  is  a  parameter 
space  for  which  the  operating  point  is  locally  stable  but  unstable  to  sufficiently 
large  perturbations. 

The  information  in  the  bifurcation  set  of  figure  11  can  be  conveyed  in  terms  of 
engine  performance.  Each  point  in  the  parameter  set  corresponds  to  a  different 
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Figure  6:  Speedline  for  fishaft  =  942s"1  (9000  RPM).  Critical  flutter  mode  is 
first  nodal  diameter,  forward  travelling  ( ND1+ ).  Critical  stall  mode  is  also  first 
nodal  diameter.  Onset  of  flutter  and  stall  instabilities  are  indicated  by  □  and 
o,  respectively. 

equilibrium  point.  In  figure  12  are  plotted  the  nondimensional  exit  pressure 
and  flow  coefficient  corresponding  to  the  stability  boundaries  of  figure  11.  In 
figure  13  are  plotted  stability  boundaries  in  terms  of  the  more  conventional 
performance  measures  of  pressure  ratio  and  inlet  Mach  number.  The  conversions 
from  nondimensional  exit  pressure  and  flow  coefficient  to  pressure  ratio  and  inlet 
Mach  number  assume  engine  dimensions  comparable  to  a  laboratory-scale  test 
facility. 
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Figure  7:  Root  locus  of  least-stable  eigenvalues  corresponding  to  speedline,  x, 
first  structural  flutter  mode,  o,  second  structural  flutter  mode. 
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Figure  9:  Mode  1  flutter  damping  by  nodal  diameter,  $  =  0.65,  =  16. 
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Figure  11:  Two  parameter  bifurcation  diagram.  Linear  stability  boundaries  as 
throttle  and  shaft  speed  are  independently  varied.  The  critical  boundaries  are 
indicated  with  heavy  lines;  —  NDO  flutter  boundary,  -  -  ND1+  flutter  boundary. 
Operation  in  the  region  to  the  right  (higher  values  of  throttle)  of  the  critical 
stability  boundaries  is  stable. 
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Figure  12:  Stability  boundaries  as  throttle  and  shaft  speed  are  independently 
varied.  The  critical  boundaries  are  indicated  with  heavy  lines;  —  NDO  flutter 
boundary,  —  ND1+  flutter  boundary.  Operation  in  the  region  to  the  right 
(higher  values  of  flow  coefficient)  of  the  critical  stability  boundaries  is  stable. 
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Figure  13:  Stability  boundaries  as  throttle  and  shaft  speed  are  independently 
varied.  Inlet  Mach  number  is  computed  assuming  Rnom  =  23.4  cm  and  sonic 
velocity  c  =  330  m/s.  The  critical  boundaries  are  indicated  with  heavy  lines;  — 
NDO  flutter  boundary,  -  -  ND1+  flutter  boundary.  Operation  in  the  region  to 
the  right  (higher  values  of  inlet  Mach  number)  of  the  critical  stability  boundaries 
is  stable. 
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6  Parametric  Sensitivity 

The  sensitivity  of  system  stability  to  parameter  variations  is  examined  by  com¬ 
puting  partial  derivatives  of  the  damping  ratios  of  several  flutter  and  stall  modes 
with  respect  to  system  parameters  (table  2).  A  point  on  the  speedline  $  =  0.55 
near  the  peak  pressure  rise  is  taken  as  the  reference.  Nine  flutter  modes  were 
examined,  zero  nodal  diameter,  and  the  first  four  nodal  diameters,  forward  and 
backward  travelling  modes.  The  first  four  stall  modes  were  examined.  The 
parametric  sensitivity  will  vary  of  course  with  the  choice  of  equilibrium  and 
nominal  parameter  values. 

The  sensitivity  of  first  nodal  diameter,  forward  travelling  flutter  damping  to 
parameter  variations  is  indicated  in  figure  14.  The  other  flutter  modes  show 
greatest  sensitivity  to  the  same  parameter  set  although  in  varying  degree  and 
sense  (positive  or  negative).  Not  surprisingly  flutter  damping  was  found  to 
be  most  strongly  influenced  by  the  structural  damping  parameters,  £  and  Ct* 
Structural  damping  is  stabilizing  to  all  of  the  flutter  modes.  Next  in  order  of 
importance  are  the  parameters  which  govern  the  natural  frequencies  and  the 
plunge  to  twist  ratios,  fea,  fcg,  e,  Qt  and  Qb-  Also  important  is  the  center 
of  pressure  offset  parameter  £cp  which  strongly  affects  the  gain  (section  2.6). 
In  contrast,  the  chord  and  blade  inertia  parameters  cr  and  D  affect  the  gain 
only  weakly.  As  discussed  in  section  5  flutter  damping  increases  with  throttle 
parameter  A.  Flutter  damping  is  shown  to  be  sensitive  to  the  orientation  and 
shape  of  both  rotor  and  stator  as  well  as  the  loss  and  deviation. 

The  sensitivity  of  first  nodal  diameter  stall  damping  to  parameter  variations 
is  indicated  in  figure  15.  The  other  stall  modes  are  sensitive  to  the  same  param¬ 
eters  and  for  some  of  the  parameters  the  sense  of  the  dependence  is  uniform.  As 
expected,  the  damping  of  the  stall  modes  is  most  sensitive  to  those  parameters 
which  determine  the  steady  characteristic  and  to  the  throttle  parameter.  For 
these  parameters,  variations  are  either  stabilizing  for  all  modes,  or  de-stabilizing 
for  all  modes.  For  other  parameters  which  strongly  influence  stall  damping,  vari¬ 
ations  may  be  either  stabilizing  or  destabilizing  depending  on  the  stall  mode. 
These  parameters  include  the  lagged  loss  parameters  rr  and  rs,  the  parameters 
which  determine  the  inertia  of  fluid  within  the  blade  passages  cs/  cos7s  and  cr, 
the  parameters  which  determine  the  inertia  of  fluid  in  the  ducts  c*iniet  and  Libs , 
and  the  parameters  which  determine  the  plunge  to  twist  ratio  and  frequency. 

6.1  Sensitivity  of  flutter  to  ratio  of  plunge  to  twist. 

Apart  from  the  structural  damping  parameters,  which  directly  affect  flutter 
damping,  the  most  influential  parameters  to  flutter  damping  are  those  which 
govern  the  plunge  to  twist  ratio  of  blade  deformation  its  natural  frequencies. 
Since  aero-loading  and  structural  damping  are  small,  the  frequencies  and  plunge 
to  twist  ratios  are  determined  primarily  by  the  conservative  form  of  the  blade 
deflection  equations. 
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Figure  14:  Parametric  sensitivity  of  ND1+  flutter  damping.  Values  plotted  are 
partial  derivatives  of  the  damping  ratio  with  respect  to  the  indicated  parameter. 


The  ratios  of  plunge  to  twist  corresponding  to  the  natural  frequencies  loi 
and  u)t  (equation  32)  are  given  by 


±\  (jeaJ^jcgWl 
Cali  Iea[w?(l  -H)-Qt] 
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Figure  15:  Parametric  sensitivity  of  ND1  stall  damping.  Values  plotted  are 
partial  derivatives  of  the  damping  ratio  with  respect  to  the  indicated  parameter. 


JL  —  (£ea  —  t.cg)DQ\ 

ca  t  ~  Ie*[u?(l-n)  -Qt)' 

respectively. 

The  sensitivity  of  the  relative  stability  of  the  flutter  modes  to  the  plunge- 
to-twist  ratio  is  illustrated  in  figure  16  for  a  small  value  of  /i  and  figure  17  for 
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a  value  of  \i  near  its  upper  limit.  Parameters  Qb  and  Qt  are  adjusted  so  that 
the  natural  frequencies  remain  equal  to  their  nominal  values.  The  highest  nodal 
diameter  patterns  are  found  to  be  the  least  stable  for  large  positive  values  of 
the  plunge-to-twist  ratio.  As  the  ratio  is  decreased,  the  critical  nodal  diameter 
pattern  decreases  to  the  zero  flutter  harmonic.  The  interblade  phase  angle  a  is 
related  to  the  flutter  harmonic  number  n  by 


a 


27r  n 


In  the  present  case  Nb  =  16,  hence  figures  16  and  17  range  in  interblade  phase 
angle  from  -180  to  180  degrees.  The  two  figures  are  qualitatively  identical 
indicating  that  \i  has  little  influence  on  the  relative  stability  of  the  flutter  modes. 


Figure  16:  Relative  stability  of  the  flutter  modes  as  the  ratio  of  the  plunge  to 
twist  is  varied  for  //  small.  The  numbers  by  each  curve  indicates  the  value  of 
q/ca ,  the  plunge  to  twist  ratio.  ax  =  0.19,  ub  =  1.81,  and  ut  =  4.60. 

It  is  shown  in  figure  18  that  structural  damping  affects  all  of  the  flutter 
modes  uniformly.  The  relative  damping  between  the  modes  is  unchanged. 

6.2  Sensitivity  of  stability  to  slope  of  the  speedline 

Stall  and  surge  instability  are  always  found  to  be  near  the  peak  of  the  total 
to  static  pressure  rise.  Gysling  [4]  developed  the  argument  that  instability 
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Figure  17:  Relative  stability  of  the  flutter  modes  as  the  ratio  of  the  plunge  to 
twist  is  varied  for  fi  large.  The  numbers  by  each  curve  indicates  the  value  of 
q/ca ,  the  plunge  to  twist  ratio.  ^/^max  =  0.93,  =  1.81,  and  u;t  =  4.60. 
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Figure  18:  Influence  of  structural  damping  on  flutter  stability.  Both  structural 
damping  parameters  are  adjusted  concurrently,  =  Ct  =  C  The  measured 
values  of  rotor-frame  frequency  are  taken  to  be  =  1.81,  and  u>*  =  4.60  and 
ft/ ft  max  ~  0.93. 
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occurs  when  mechanical  energy  is  introduced  into  the  flow  faster  than  it  can  be 
convected  away.  The  mechanical  energy  introduced  into  the  flow  is  proportional 
to  the  total  to  total  pressure  rise, 


assuming  an  exit  flow  angle  of  zero.  Conservation  of  mass  requires  that 

$in  =  *out  = 

The  kinetic  energy  of  the  flow  can  be  convected  downstream,  hence  the  energy 
balance  is  unstable  to  perturbations  if 

OT«>a(y), 


or 


d9  tt 


>$. 


The  instability  criterion  may  be  re-written 


d9 ts 


>0, 


since  the  total  to  total  pressure  rise  is  related  to  the  total  to  static  pressure  rise 
by 

$2 

*  tt  =  *ts  +  *y* 

High-incidence  flutter  instability  is  also  typically  found  near  the  peak  of 
the  compressor  characteristic.  Figures  14  and  15  indicate  that  variations  in 
parameters  which  affect  the  steady  compressor  characteristic  have  the  same 
effect  on  the  stability  of  flutter  as  stall.  The  degree  of  flutter  sensitivity  is  far 
less  however. 

Damping  of  the  flutter  modes  is  found  to  be  somewhat  dependent  on  the 
parameters  governing  loss  and  deviation.  In  figure  19  are  plotted  values  of  the 
flutter  damping  ratios  for  zero  loss  and  deviation,  nominal  loss  and  deviation, 
and  loss  and  deviation  fifty  percent  greater  than  nominal.  All  other  parameters 
are  fixed  at  their  nominal  values.  The  relative  damping  between  flutter  modes  is 
clearly  shown  to  be  dependent  on  viscous  effects;  the  least-damped  flutter  mode 
shifts  from  third  nodal  diameter  forward-travelling  to  the  zero  nodal  diameter 
mode  as  losses  and  deviation  increase.  Flutter  is  more  damped  for  the  nominal 
loss  and  deviation  than  for  either  of  the  other  cases.  A  significant  implication 
is  that  inviscid  computational  fluid  dynamics  codes  cannot  accurately  predict 
the  onset  of  flutter  if  there  is  no  additional  empirical  accounting  for  loss  and 
deviation.  Such  codes  will  be  most  accurate  at  predicting  flutter  damping  for 
conditions  of  small  loss  and  deviation. 
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Figure  19:  Sensitivity  of  flutter  damping  to  loss  and  deviation.  Damping  ratios 
are  plotted  corresponding  to  scaling  of  the  loss  and  deviation  by  a  factors  of  0 
o,  l  x,  and  1.5  o.  The  flow  coefficient  for  these  results  is  $  =  0.55.  All  other 
parameters  take  their  nominal  values. 


If  the  losses  are  increased,  the  compressor  peak  will  move  right  (toward 
higher  values  of  flow  coefficient)  and  the  compressor  will  be  less  stable.  Clearly 
then,  increases  in  the  loss  coefficients  Lr i,  Lr 2,  Lr 3,  Ls  1,  Ls 2,  and  Ls 3  destabilizes 
the  system  because  the  incidence  angles  are  nearly  always  positive.  Because  the 
incidence  angles  are  less  than  unity,  Lr 3  has  greater  influence  than  Ir2  which 
in  turn  has  greater  influence  than  Lr  1 .  A  similar  relation  holds  for  the  stator 
loss  coefficients.  Increases  in  the  zero-incidence  angles  of  the  blades  /?zr  and  /?zs 
decreases  the  incidence  angles  (figure  4)  which  decreases  the  loss  and  improves 
compressor  stability. 

Increased  turning  of  the  flow  results  in  greater  pressure  rise.  If  the  pressure 
rise  is  increased  the  compressor  peak  will  move  towards  lower  values  of  the  flow 
coefficient  increasing  system  stability.  This  can  be  accomplished  by  decreasing 
the  trailing-edge  metal  angle  f3r  and  the  flow  deviation  parameters  Ai  and  A2. 

One  would  suppose  then  that  blades  would  be  designed  with  large  values  of 
the  zero-incidence  angle  and  small  values  of  the  trailing-edge  metal  angle.  This 
underscores  the  limitations  of  the  model.  In  this  model  the  loss  and  deviation 
coefficients  enter  as  parameters  which  are  to  be  identified  empirically.  The 
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implicit  dependence  on  blade  geometry  is  hidden.  To  a  great  extent,  these 
empirical  parameters  are  determined  from  the  blade  shape  and  it  is  known  that 
extreme  curvature  of  the  blades  can  greatly  increase  the  loss. 

6.3  Sensitivity  of  stall  to  blade  deflection  parameters 

Steady  blade  deformation  affects  the  characteristic  in  two  ways.  An  increase 
in  the  blade  twist  increases  the  incidence  angle,  and  hence  the  loss,  which  is 
de-stabilizing.  Positive  blade  twist  also  increases  flow  turning  which  is  a  stabi¬ 
lizing  effect.  For  this  set  of  parameters,  the  steady  blade  deflection  is  positive. 
Hence,  all  parameter  variations  which  increase  the  steady  twist  will  tend  to  af¬ 
fect  damping  in  the  same  manner.  Increases  to  the  torsional  stiffness  Qt>  the 
center  of  pressure  offset  £cp,  and  the  rotational  inertia  e  all  have  the  effect  of 
decreasing  the  steady  twist.  For  stall  modes  1,3,  and  4,  these  variations  increase 
the  damping  but  for  mode  2,  the  damping  is  decreased  (figure  20).  Increases  to 
the  elastic  axis  offset  £ea  and  the  center  of  pressure  offset  £cp  increase  the  gain 
of  the  aeroload  and  are  found  to  de-stabilize  modes  1,3,  and  4,  and  stabilize 
mode  2. 

The  anomalous  behavior  of  mode  2  can  only  be  explained  by  unsteady  in¬ 
teraction  between  blade  deflection  and  flow  perturbations.  This  is  supported 
by  the  observation  that  the  blade  damping  parameters  have  very  little  influence 
on  the  stability  of  modes  1,3,  and  4  while  the  twist  damping  parameter  Ct  does 
have  some  influence  on  the  stability  of  mode  2. 

The  parameters  with  greatest  influence  are  £cp,  £Cg>  and  e.  The  rotor  chord 
cr  also  is  found  to  have  a  strong  influence  on  stall  damping,  but  the  dependence 
is  more  complicated;  it  contributes  to  the  gain  of  the  aeroload  and  the  inertia 
of  the  fluid  in  the  blade  passage. 

6.4  Blade  stiffness 

The  stiffness  parameters  <2b  and  Qt  are  not  among  the  most  influential  with 
respect  to  small  deviations  about  the  nominal  parameter  set.  As  one  considers 
very  large  values  however,  it  can  be  shown  that  the  flutter  modes  are  progres¬ 
sively  stabilized.  Finally,  the  blades  are  effectively  rigid  and  the  only  instability 
modes  are  stall  and  surge. 

In  figure  21  are  plotted  the  damping  ratios  of  the  flutter  modes  by  nodal 
diameter  pattern  for  progressively  higher  values  of  the  stiffness  parameters.  In¬ 
creasing  stiffness  is  shown  to  stabilize  some  flutter  modes  and  to  destabilize 
others  through  reduction  of  the  aerodamping.  Eventually,  the  aerodamping  is 
negligible  in  comparison  to  the  structural  damping  and  all  modes  are  equally 
damped. 
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Figure  20:  Parametric  sensitivity  of  stall  damping  to  parameters  governing  blade 
deflection.  Values  plotted  are  partial  derivatives  of  the  damping  ratio  with 
respect  to  the  indicated  parameter. 


6.5  Unsteady  loss 

The  effect  of  modeling  unsteady  losses  with  first  order  lag  equations  on  stall 
stability  has  been  reviewed  by  Longley  [4].  In  this  section,  those  observations 
are  verified  for  the  present  model  and  the  interpretation  presented  by  Longley 
is  briefly  re-capitulated. 

The  stability  of  the  rotating  stall  modes  is  clearly  affected  by  blade  deforma- 
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Figure  21:  Relative  stability  of  the  flutter  modes  as  the  blade  stiffness  is  varied. 
The  numbers  by  each  curve  refer  to  the  value  of  Qb\  2.78  is  its  nominal  value. 
The  ratio  Qt/Qb  is  held  constant. 


tion  and  by  the  inclusion  of  the  inter-bladerow  space.  To  isolate  the  influence 
of  the  time  lags,  the  blade  stiffness  parameters  are  greatly  increased  and  the 
length  of  the  inter-bladerow  space  is  greatly  decreased. 

For  zero  lag,  damping  of  all  stall  modes  is  equal  at  the  flow  coefficient  corre¬ 
sponding  to  the  compressor  peak  as  shown  in  figure  23.  That  value  of  damping 
is  positive  however  and  instability  occurs  in  the  first  stall  mode  just  beyond 
peak.  This  discrepancy  is  probably  due  to  differences  in  the  present  formula¬ 
tion  of  the  rotor  energy  equation  and  the  unsteady  Bernoulli  equation  which  is 
conventional  for  modeling  stall. 

Consistent  with  the  results  presented  by  Longley,  nonzero  lag  is  found  to  be 
stabilizing  for  low  flow  rates,  figure  25.  The  amount  of  stabilization  is  greater 
for  the  higher  modes.  At  the  compressor  peak  however,  damping  ratios  for 
the  first  and  second  stall  modes  are  actually  lower  than  those  calculated  with 
steady  losses.  Comparing  figures  22  and  24  it  is  clear  that  lagging  of  the  losses 
results  in  higher  frequencies  than  quasi-steady  losses  which  is  also  consistent 
with  Longley. 

Longley  has  summarized  earlier  work  demonstrating  that  in  simplified  mod¬ 
els  such  as  this  one,  lagged  losses  and  deviation  increase  the. effective  inertia  of 
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Figure  22:  Root  locus  of  stall  eigenvalues  for  steady  loss,  short  inter-bladerow 
space  and  effectively  rigid  blades.  —  1st  mode,  —  2nd  mode,  -  •  3rd  mode,  •  •  • 
4th  mode.  All  other  parameters  take  their  nominal  values. 


fluid  proportionately  to  the  stall  mode  number  within  the  blade  passages.  The 
argument  was  derived  for  compressors  and  not  for  fans  and  does  not  incorporate 
inter-bladerow  flow  domain.  Alternatively,  one  may  think  of  unsteady  losses  as 
decreasing  the  instantaneous  slope  of  the  compressor  characteristic. 

Blade  deformation  and  the  presence  of  the  inter-bladerow  gap  dramatically 
alter  the  effect  of  unsteady  losses  as  shown  in  figure  26.  In  contrast  to  the 
short-gap,  rigid-blade  case,  lagged  losses  are  shown  to  de-st abilize  the  third  and 
fourth  mode  and  stabilize  the  first  and  second  modes.  Further,  the  root  locus 
of  the  second  mode  is  greatly  affected.  In  both  cases  however,  lagging  of  the 
losses  results  in  an  increase  in  the  stall  frequency.  This  complicated  interplay 
of  parametric  sensitivity  warrants  further  attention.  In  particular,  it  must  be 
verified  that  the  Pade  approximant  is  not  introducing  non-physical  effects. 

The  effect  of  loss  time  lags  on  flutter  stability  is  illustrated  in  figure  27.  Loss 
time  lags  are  shown  to  have  a  stabilizing  influence  on  the  forward-travelling 
waves  and  a  destabilizing  influence  on  the  backward-travelling  waves.  The  flut¬ 
ter  modes  are  far  less  sensitive  than  are  the  stall  modes,  however  the  flutter 
modes  are  also  lightly  damped  and  such  differences  may  be  very  significant. 
Certainly,  this  result  may  depend  strongly  on  the  choice  of  other  parameters. 
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Figure  23:  Damping  ratios  of  stall  for  steady  loss,  short  inter-bladerow  space 
and  effectively  rigid  blades.  —  1st  mode,  —  2nd  mode,  -  •  3rd  mode,  •  •  •  4th 
mode.  All  other  parameters  take  their  nominal  values. 

6.6  Installation  effects 

The  inlet  geometry  parameter  and  the  length  of  the  inter-bladerow  space  can 
effect  the  stability  of  stall  and  flutter  modes.  As  discussed  in  section  2.1,  inlet 
geometry  clearly  effects  the  wave  numbers  of  the  flow  perturbations.  Another 
effect  of  inlet  geometry  is  that  the  flow  velocity  must  be  straightened  as  it 
enters  the  duct;  a  longer  inlet  tends  to  straighten  the  flow  better  upstream 
of  the  rotor  face.  The  present  model  assumes  nominally  two-dimensional  flow 
and  cannot  capture  this  effect.  Within  an  inter-bladerow  space  it  has  been 
shown  in  section  2.7  that  the  perturbation  to  the  flow  field  will  have  potential 
components  (the  homogenous  solution)  which  decay  exponentially  from  each 
end  of  the  domain  and  a  vortical  component  which  is  convected  downstream. 
A  long  inter-bladerow  space  allows  the  potential  perturbation  to  significantly 
decay,  partially  decoupling  the  bladerows.  The  effect  of  inter-bladerow  gaps  on 
stall  instability  is  discussed  in  the  review  of  Longley  [4].  In  such  a  case  one 
must  consider  the  total  to  static  pressure  rise  of  each  bladerow.  Stall  instability 
is  then  found  to  occur  when  a  single  bladerow  characteristic  has  zero  slope  and 
the  others  all  have  negative  slope. 

The  effect  of  the  inlet  geometry  parameter  on  stall  and  flutter  stability  is 
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Figure  24:  Root  locus  of  stall  eigenvalues  for  lagged  loss,  short  inter-bladerow 
space  and  effectively  rigid  blades.  —  1st  mode,  —  2nd  mode,  -  •  3rd  mode,  •  *  • 
4th  mode.  All  other  parameters  take  their  nominal  values. 

illustrated  in  figures  28  and  29.  The  higher-order  stall  modes  are  found  to 
decrease  in  frequency  and  are  significantly  damped  as  the  inlet  parameter  ainiet 
increases  from  0,  an  infinite  inlet  duct,  to  2,  a  very  short  inlet  duct.  The  first 
stall  mode,  which  is  the  critical  mode  for  instability,  decreases  in  frequency 
but  its  damping  ratio  is  nearly  unchanged.  The  inlet  parameter  is  shown  to 
have  no  influence  on  flutter  modes  corresponding  to  nodal  diameters  0,  +8, 
and  —8.  The  other  forward  travelling  modes  are  found  increase  in  damping  as 
ainiet  increases  (decreasing  duct  length);  the  other  backward  travelling  waves 
are  found  to  decrease  in  damping. 

The  dependence  of  flutter  instability  on  the  length  of  the  inter-bladerow 
space  is  illustrated  in  figure  30.  As  in  the  case  of  inlet  geometry  variations,  the 
inlet  parameter  is  shown  to  have  no  influence  on  flutter  modes  corresponding 
to  nodal  diameters  0,  +8,  and  -8.  In  contrast  to  the  inlet  duct  however,  long 
lengths  of  the  inter-bladerow  space  are  found  to  be  de-stabilizing  to  both  forward 
and  backward-travelling  modes. 

The  dependence  of  stall  instability  on  the  length  of  the  inter-bladerow  space 
is  much  more  complicated.  As  shown  in  figure  31,  there  appear  to  be  multiple 
root  loci  for  each  of  the  stall  modes.  This  may  be  indicative  that  the  nature 
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Figure  25:  Damping  ratios  of  stall  for  lagged  loss,  short  inter-bladerow  space 
and  effectively  rigid  blades.  —  1st  mode,  —  2nd  mode,  -  •  3rd  mode,  •  •  •  4th 
mode.  All  other  parameters  take  their  nominal  values. 


of  the  stall  mode  changes  as  the  bladerows  become  effectively  de-coupled  or  it 
could  simply  be  an  artifact  of  the  Pade  approximant  used  in  the  modeling  of  the 
perturbation  flow  field  (section  2.7).  For  large  values  of  Zibs,  the  stall  modes  are 
found  to  have  much  higher  frequencies  and  less  damping.  Though  not  shown 
here,  the  root  loci  of  the  flutter  modes  are  found  to  be  continuous  over  this 
range  of  Libs- 

6.7  Inlet  total  pressure 

The  nondimensional  inlet  total  pressure  will  vary  with  operating  conditions.  In 
particular,  it  increases  with  both  airspeed  and  temperature.  Flutter  is  sensitive 
to  inlet  total  pressure  but  stall  is  not.  The  dependence  of  flutter  damping  on 
inlet  total  pressure  is  illustrated  in  figure  32.  Increased  inlet  total  pressure 
is  shown  to  destabilize  the  forward-traveling  modes  and  the  first  backward¬ 
traveling  mode,  the  other  backward-traveling  modes  are  stabilized,  and  there  is 
no  influence  on  the  eighth  nodal  diameter  mode. 
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Figure  26:  Root  locus  of  stall  eigenvalues  as  flow  coefficient  varies  from  0.4  to 
0.8  for  lagged  loss  x  and  for  steady  loss  o.  All  other  parameters  take  their 
nominal  values. 
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Figure  27:  Sensitivity  of  flutter  damping  to  unsteady  losses;  lagged  loss  x, 
steady  loss  o.  The  flow  coefficient  is  given  by  $  =  0.55.  All  other  parameters 
take  their  nominal  values. 
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Figure  29:  Effect  of  c^niet  on  damping  ratios  of  flutter  modes.  ainjet  takes  the 
values  0  o,  0.5  □,  1  x,  and  2  o.  The  flow  coefficient  is  given  by  $  =  0.55.  All 
other  parameters  take  their  nominal  values. 
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Figure  30:  Dependence  of  flutter  damping  ratios  on  length  of  the  inter-bladerow 
space,  Libs-  ^ibs  takes  the  values  2  x  10“3  o,  0.1  □,  0.2  x,  and  0.5  o.  The  flow 
coefficient  is  given  by  $  =  0.55.  All  other  parameters  take  their  nominal  values. 
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Figure  31:  Root  locus  of  stall  eigenvalues  as  length  of  the  inter-bladerow  space 
is  varied.  Libs  varies  from  2  x  lO”3,  o,  to  0.5,  □.  The  flow  coefficient  is  given  by 
$  =  0.55.  Blades  are  effectively  rigid,  all  other  parameters  take  their  nominal 
values. 
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Figure  32:  Dependence  of  flutter  damping  ratios  on  inlet  total  pressure,  Patm- 
Patm  takes  the  values  1.5  o,  2.0  x,  2.5  □,  and  3.0  o.  The  flow  coefficient  is  given 
by  $  =  0.55.  All  other  parameters  take  their  nominal  values. 
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7  System  identification  and  control 

Experimentally,  system  stability  may  be  evaluated  from  its  forced  response.  In 
this  section,  the  role  of  actuation  in  system  identification  is  investigated.  Mass 
flow  actuation  will  be  used  to  drive  the  system,  system  output  which  will  be 
examined  include  pressure  at  the  rotor  trailing  edge,  blade  deflection  in  the  rotor 
frame  (such  as  might  be  measured  with  a  strain  gage),  and  blade  deflection  in 
the  stator  frame. 

The  DAE  system  is  of  the  form 

A(z,i,0  =  0, 

where  £  represents  the  modal  amplitudes  of  the  bleed  mass  flow.  For  a  par¬ 
ticular  choice  of  the  actuation  £,  the  system  transfer  functions  can  be  found 
by  first  linearizing  the  system  and  assuming  a  time  dependence  of  the  form 
x  =  z0exp(s*), 

[sAx  +  Ax]  x  +  Atf  =  0. 

This  is  a  linear  algebra  problem  of  the  form  Ax  =  b  where  A  =  sAx  +  Ax  and 
b  =z  .  The  matrix  A  may  be  singular  in  which  case  a  minimum  least- 

squares  solution  may  be  found  using  the  singular  value  decomposition.  If  the 
actuation  vector  £  is  chosen  to  have  unit  magnitude,  then  the  transfer  functions 
are  given  by  the  solution  x. 

The  transfer  functions  inherit  the  linear  independence  of  the  different  az¬ 
imuthal  harmonics.  For  example,  consider  the  static  pressure  downstream  P2a 
from  the  rotor  as  the  output.  The  modal  expansions  of  the  bleed  mass  flow  rate 
and  pressure  P^a  are  then 

N 

Z(0,  t)  =  £ori  +  22  6r/  cos(k0)  -  (kim  sin (kd) 
k-l 
'  N 

F*2a(0>O  =  P2a,0rl  4*  ^  ^  P2a,krl  COS (k9)  —  P2a,kim  sin(fcfl). 

k~l 

In  the  case  of  modal  actuation,  only  a  single  pair  of  amplitudes  tkrh  €kim 
is  nonzero.  It  can  easily  be  shown  that  the  output  responds  only  in  that  same 
mode,  p2a}jrl  =  P2a}jim  =  0,  for  all  j  ^  k.  Further,  there  is  only  one  indepen¬ 
dent  transfer  complex  function  Gk  between  the  inputs  and  outputs: 

P2a,krl  *4*  i-^2 a,kim  ri  /  \ 

- 7 - —7 -  ”  Gk(S )• 

S krl  “r  Isfcim 

Point-to-point  transfer  functions  refer  to  the  response  of  a  particular  sensor 
to  a  particular  actuator.  In  this  case,  a  single  bleed  valve  located  at  6  —  Of,  is 
actuated, 

Z(e,t)  =  tPe‘i6(e-eb). 
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The  delta  function  may  be  expanded  in  a  Fourier  series, 
tpe9iS(O-0b)  =  tpe9i 

Superposition  applies  and  the  response  of  a  particular  pressure  sensor  at  9  —  9p 
is  found  to  be 

p(9p ,  t)  =  c 

in  terms  of  the  modal  transfer  functions. 

Of  particular  interest  is  the  transfer  function  between  blade  deflection  and 
a  travelling  wave  generated  by  bleed  actuation.  The  flutter  and  stall  modes  are 
both  travelling- waves.  Travelling- wave  excitation  would  take  the  form, 

£(0,0  =&  cos (k9  +  ut  +  <p). 

In  the  present  model  formulation,  the  rotor  has  angular  velocity  -1.  Forward¬ 
travelling  waves  travel  in  the  same  direction  as  the  rotor,  hence,  k  >  0  corre¬ 
sponds  to  forward-travelling  waves  and  k  <  0  corresponds  to  backward  travelling 
waves.  A  transfer  function  Gfc(u>)  can  be  found  such  that  the  blade  deflection 
(measured  in  the  stator  frame)  will  be  given  by 

a(9,t)  =  Gk(u)£k  cos (k9  +  art  4-  <j>). 

If  xp  is  taken  to  be  the  azimuthal  coordinate  fixed  to  the  rotor  frame,  the  coor¬ 
dinate  transformation  will  be  given  by  9  =  xp  -  t.  The  blade  deflection  in  the 
rotor  frame,  which  could  be  measured  with  a  strain  gage,  is  then  given  by 

a v{xp,t)  =  cos [kip  +  (w  -  k)t  +  <p]. 

The  transfer  function  between  any  output  and  any  input  contains  stability 
information  about  all  of  the  system  eigenmodes.  In  practice  however,  it  is  im¬ 
portant  to  choose  actuation  that  acts  strongly  upon  the  eigenmode  of  interest. 
Otherwise  the  information  may  be  lost  in  experimental  noise.  It  is  also  impor¬ 
tant  to  select  sensors  for  which  the  direct  effect  of  actuation  is  comparable  to  or 
smaller  than  the  response  due  to  excitation  of  the  system  mode.  This  is  evident 
in  the  present  model.  Blade  deflection  is  clearly  a  good  measure  of  the  flutter 
instability.  Suppose  blade  deflection  and  mass  flow  actuation  are  related  by  a 
transfer  function  G, 

a  =  G(s)£. 

If  G  were  approximated  by  a  ratio  of  polynomials,  the  roots  of  the  numerator, 
the  zeroes ,  would  be  well-separated  from  the  roots  of  the  denominator,  the  poles. 
This  implies  that  if  the  excitation  is  very  small,  the  response  might  still  be  large 


r  i  N 

77^  +  -  Gk  cos (k9b  -  k9p) 
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in  the  neighborhood  of  the  poles.  The  poles  are  identical  to  the  eigenvalues  of 
the  unforced  system.  The  experimentally  measured  transfer  function  can  then 
be  used  to  estimate  system  eigenvalues.  Alternatively,  consider  how  system 
stability  might  be  interpreted  from  the  transfer  function  between  trailing-edge 
pressure  and  mass  flow  actuation.  Blade  deflection  alone  will  influence  pressure 
to  some  degree;  mass  flow  actuation  alone  will  have  a  much  greater  impact,  thus 

p  =  Aa  + 

where  B  >  A  over  the  frequency  range  of  interest.  The  measured  transfer 
function  between  pressure  and  mass  flow  actuation  will  then  be 

p  =  (AG  +  B){. 

Representing  G  as  a  ratio  of  its  numerator  and  denominator  G  =  N/D> 

AN  +  BD, 

p' — — «• 

Since  B  A,  the  transfer  function  will  have  zeros  very  near  the  poles  associated 
with  flutter.  Coupled  with  noise,  this  greatly  complicates  the  task  of  estimating 
the  poles  from  the  experimentally  obtained  transfer  functions.  As  an  example, 
transfer  functions  between  blade  deflection  and  mass  flow  and  between  trailing- 
edge  pressure  and  mass  flow  for  second  nodal  diameter,  forward-travelling  flutter 
are  illustrated  in  figures  33  and  34.  The  blade  deflection  transfer  function  indi¬ 
cates  the  classic  response  of  a  forced  oscillator.  In  contrast,  the  “kink”  in  the 
pressure  transfer  function  is  consistent  with  pole-zero  cancellation.  Rational 
polynomials  have  been  fit  to  the  curves  to  illustrate  how  one  might  estimate 
system  stability  information  from  experimentally-measured  transfer  functions. 

From  the  model,  the  poles  and  zeros  may  be  computed  directly  for  either 
source  of  outputs.  The  poles,  of  course,  won't  change.  Returning  again  to  the 
linearized  DAE  system, 

[sAi  +  Ax]  x  +  Ag  =  0, 

the  poles  correspond  to  values  of  s  which  satisfy  the  linearized  system  when 
£  =  0.  Likewise,  the  zeros  correspond  to  values  of  s  which  satisfy  the  linearized 
system  when  the  output  y  =  Ax  —  0.  Each  case  can  be  solved  as  an  eigenvalue 
problem,  the  poles  are  given  by  the  eigenvalues  of  the  matrix 

[sAx  +  Ar] , 

and  the  zeroes  are  given  by  the  eigenvalues  of  the  matrix 

'  (sAx  4-  Ax)  Ac  ' 

A  0  _  * 

The  actual  values  for  the  poles  estimated  in  figures  33  and  34  are  0.00132±3.80i; 
the  actual  values  for  the  zeroes  estimated  in  figure  34  are  —0.0947  ±  3.83i. 
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Figure  33:  Transfer  function  between  blade  deflection  and  mass  flow  actuation, 
2nd  nodal  diameter,  forward-travelling  flutter  mode.  —  indicates  a  fit  using  2 
poles  and  2  zeroes.  Estimated  eigenvalues:  0.00125  ±  3.80i;  estimated  zeroes: 
7.10  and  -0.903. 
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Figure  34:  Transfer  function  between  trailing-edge  pressure  and  mass  flow  ac¬ 
tuation,  2nd  nodal  diameter,  forward-travelling  flutter  mode.  —  indicates  a  fit 
using  2  poles  and  4  zeroes.  Estimated  eigenvalues:  0.00120  ±  3.80i;  estimated 
zeroes:  0.0110  ±  3. 79i,  23.5,  and  2.35. 
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8  Passive  control 

Based  on  the  observation  that  flow  perturbations  and  blade  deflection  are  cou¬ 
pled,  it  is  plausible  that  blade  deflection  can  be  treated  indirectly  by  adding 
damping  to  the  flow  perturbations.  In  this  section,  Helmholtz  resonators  are 
examined  with  regard  to  passive  control  of  flutter.  Consistent  with  the  assumed 
form  of  mass  flow  actuation,  it  is  desirable  that  passive  dampers  act  on  all  flutter 
modes.  This  can  be  accomplished  by  distributing  Helmholtz  resonators  finely 
and  evenly  around  the  casing  as  in  figure  35. 


Figure  35:  Schematic  of  Helmholtz  resonators  for  passive  control,  1/4  annulus. 

Considering  a  single  resonator,  the  conservation  of  mass  within  the  plenum 
may  be  expressed 

cz 

where  pi  is  the  pressure  within  the  plenum,  m,-  is  the  mass  flow  into  the  plenum, 
V  is  the  volume  of  the  plenum,  and  c  is  the  sonic  velocity.  The  conservation  of 
momentum  within  the  neck  of  the  resonator  may  be  expressed 

=  “  Pi^A  ~  K™ij 

where  p(0i,t)  is  the  pressure  directly  outside  the  plenum,  /  is  the  length  of  the 
neck,  A  is  the  cross-sectional  area  of  the  neck,  and  k  is  a  viscous  loss  constant.  In 
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the  limit  of  a  continuum  of  small  oscillators,  these  equations  may  be  combined 
and  nondimensionalized  to  obtain  the  oscillator  equation 

£  +  2ChQhzi*  +  Qh£  =  PnPzb-  (33) 

The  damping  ratio  is  given  by  £h  and  Qh  is  the  ratio  of  the  Helmholtz  resonator 
frequency  to  the  shaft  frequency;  the  gain  /?h  =  A/(nRnoml ).  The  DAE  system 
(30)  is  augmented  by  equation  (33)  and  the  mass  flow  actuation  parameters  now 
become  state  variables. 


9  Mistuning 

It  has  been  noted  that  the  linearized  governing  equations  de-couple  by  the 
Fourier  harmonics.  This  is  a  consequence  of  azimuthal  symmetry.  If  that  sym¬ 
metry  is  broken,  there  will  be  some  coupling  of  the  Fourier  harmonics  in  each 
of  the  stall  and  flutter  mode.  One  might  reasonably  hypothesize  that  this  cou¬ 
pling  would  lead  to  increased  damping  in  the  least  damped  modes  and  decreased 
damping  in  the  most  damped  modes.  A  negative  aspect  of  mistuning  is  that  any 
disturbance  will  excite  all  of  the  modes  to  some  degree.  The  forced-response  of 
the  system  will  then  be  larger. 

Mistuning  has  been  incorporated  into  the  reduced-order  model  by  allowing 
for  variability  in  blade  stiffness  and  by  allowing  for  variability  in  the  stagger 
angle.  First,  the  equations  of  motion  were  re-derived  in  the  rotor  frame  of 
reference.  Secondly,  the  stagger  angle  7r  and  the  blade  stiffness  parameters  Qb 
and  Qt  were  expanded  as  Fourier  series  in  the  rotor-frame  azimuthal  coordinate 
ip.  The  Galerkin  projection  over  the  azimuthal  angle  becomes  excessively  tedious 
to  accomplish  using  symbolic  algebra.  Collocation  (as  in  the  stall  model  of 
Mansoux  et.al  [7])  has  been  used  thus  far,  but  results  are  not  satisfactory.  In 
the  future,  numerical  integration  over  ip  will  be  pursued. 

The  present  approach  to  mistuning  is  novel  in  that  it  considers  mistun¬ 
ing  of  the  complete  aeromechanical  system.  One  consequence  of  mistuning  is 
that  a  steady  solution  exists  only  in  the  mistuned  frame.  If  for  example  the 
blade  stiffness  is  mistuned,  the  steady  deflection  will  be  a  function  of  the  angle 
ip.  In  the  stator  frame,  the  observed  blade  deflection  will  be  periodic  over  a 
blade  rotation.  Another  consequence  is  that  the  blades  will  react  differently 
to  the  aeroload  and,  conversely,  the  flow  perturbations  will  be  influenced  by 
the  variable  blade  defection.  All  of  the  consequences  are  accounted  for  in  the 
reduced-order  model. 


10  Nonlinear  interaction  of  stall  and  flutter 

In  addition  to  local  stability,  one  must  also  consider  stability  with  respect  to 
finite  amplitude  perturbations.  Below  in  figure  36  are  plotted  time  histories 
of  stall  and  flutter,  respectively,  obtained  by  numerically  integrating  the  DAE 
system.  The  instabilities  occur  as  Hopf  bifurcations  and  also  break  azimuthal 
symmetry.  The  amplitude  saturation  is  due  to  nonlinearity.  As  revealed  in  these 
plots,  flutter  and  stall  are  typically  coupled. 

Rotating  stall  typically  occurs  as  a  hard  loss  of  stability  or  subcritical  Hopf 
bifurcation.  There  exist  a  range  of  parameters  near  the  stability  boundary  for 
which  the  system  is  locally  stable  but  unstable  to  finite  perturbations.  This  has 
been  observed  in  practice  and  in  experiments.  Supposing  the  engine  is  tested 
along  a  speedline,  the  system  enters  stall  near  the  compressor  peak  but  only 
recovers  from  stall  when  the  operating  point  is  returned  to  much  lower  pressure 
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Figure  36:  (a)  Time  history  of  stall  oscillations,  (b)  Time  history  of  flutter 
oscillations. 


ratios.  Thus  far,  flutter  instability  has  been  found  to  be  a  soft  loss  of  stability 
or  supercritical  Hopf  bifurcation.  A  numerically  generated  bifurcation  diagram 
is  plotted  in  figure  37  along  with  the  locus  of  fixed  points  (speedline).  The 
throttle  parameter  A  is  slowly  varied  in  time;  in  these  simulations  a  small  noise 
perturbation  has  been  added  to  the  first  harmonic  of  the  inlet  flow  coefficient. 
Initially,  the  throttle  parameter  is  gradually  decreased  from  a  suitably  large 
value  at  which  the  system  is  stable.  As  A  is  decreased,  the  equilibrium  point 
moves  up  the  speed  line  and  the  system  loses  stability  through  flutter  at  a 
flow  coefficient  of  $f  =  0.564.  The  apparent  gap  between  flutter  instability 
and  saturation  of  the  flutter  limit  cycle  is  due  to  the  fact  that  the  growth 
rate  is  initially  so  small  that  the  variation  of  A  does  not  appear  quasi-steady. 
At  $  «  0.525  the  trajectory  leaves  the  flutter  limit  and  settles  onto  the  stall 
limit  cycle.  The  throttle  parameter  is  then  gradually  increased  back  to  its 
initial  value.  Co-existent  flutter  and  stall  limit  cycles  are  demonstrated  between 
<$>  —  0.52  and  $  =  0.55.  The  turning  point  of  the  stall  limit  cycle  is  apparently 
at  <3>  «  0.55,  well  above  the  point  of  linear  instability  at  3>s  =  0.536.  Below 
this  turning  point,  the  stall  limit  cycle  does  not  exist  and  the  trajectory  quickly 
drops  back  onto  the  equilibrium  point  before  gradually  cycling  back  up  to  the 
flutter  limit  cycle.  This  plot  demonstrates  hysteresis  between  flutter  and  stall 
limit  cycles.  If  the  flutter  instability  were  not  present,  hysteresis  between  the 
stall  limit  cycle  and  the  equilibrium  curve  would  exist  in  the  range  between  the 
point  of  stall  instability  and  the  turning  point  of  the  stall  limit  cycle. 

In  future  work,  normal  forms  of  the  dynamical  system  will  be  interpreted 
from  simulations  following  the  example  of  Anderson  el  al  [1].  This  will  allow 
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Figure  37:  Numerical  bifuration  diagram  demonstrating  hysteresis  between  stall 
and  flutter  limit  cycles  due  to  subcritical  Hopf  bifurcation,  x  indicates  point  of 
flutter  instability,  o  indicates  point  of  stall  instability. 


an  analytical  approximation  of  the  unstable  branch  of  the  stall  limit  cycle  and 
analysis  of  codimension-2  bifurcations. 

11  Conclusion 

A  reduced-order  model  representing  a  unified  view  of  the  stability  of  flexibly 
bladed  compression  systems  has  been  described.  The  model  captures  the  es¬ 
sential  nonlinear  coupling  between  surge,  rotating  stall  and  flutter,  and  pro¬ 
vides  a  vehicle  for  the  systematic  investigation  of  their  interaction.  The  ap¬ 
proach  to  model  construction  described  above  provides  a  tractable  qualitative 
yet  physically-based  dynamical  description  of  the  interaction  of  performance- 
limiting  compression  system  instability  phenomena.  Analysis  of  the  linear  sta¬ 
bility  boundaries  is  consistent  with  operability  limits  observed  in  practice  for 
turbomachinery.  The  post  stall  and  flutter  behavior  simulated  is  also  consistent 
with  experimentally  observed  behaviors. 
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